SPE 115929 Modeling Leakage through Faults of CO2 Stored in an Aquifer

For secure storage of CO2 within geologic formations, the integrity of caps – overlying strata that are impervious to CO2 – is an important factor. Geologic structures, notably faults and the damage zones surrounding them may provide a conduit for CO2 to escape through a cap. If the fault encounters shallower permeable formations, the CO2 rising along the fault can enter them. This lateral migration would attenuate the rate at which CO2 enters sensitive formations such as aquifers used for drinking water. Thus, CO2 leakage along faults will have three behaviors: upward migration from the storage formation along a fault, lateral movement from the fault into permeable layers, and a continued but attenuated CO2 flux along the fault above the layers. Here we develop a quasi-1D single-phase flow model for these three behaviors. The model is highly simplified and intended to be suitable for inclusion in a certification framework for geologic storage projects. The model accounts for flow from the fault into a permeable formation using a leakoff coefficient. The coefficient can vary spatially and depends on the geometry and petrophysical properties of the formation. We apply a commercial simulator to verify the quasi-1D model. A series of examples illustrates the controlling mechanisms for leakage rate from the reservoir and its attenuation by flux into shallower layers. Nonlinearities arise even in this simple model. For example, leakage flux and the degree of attenuation vary nonlinearly with the permeability of the fault and the permeability of the shallower layer(s) intersected by the fault. Layers nearest the CO2 storage formation produce the most attenuation. But the percentage of CO2 entering overlying formations from the fault varies linearly with the ratio of fault permeability to leakoff coefficient. A simple estimate of the leak-off coefficient compares favorably with 2D, full-physics simulations. If the permeable layer is dipping, CO2 enters it asymmetrically and estimating the leakoff coefficient is less straightforward. The difference arises because of preferential flow within the layer (CO2 in the upper part, water below). Introduction If society elects to reduce anthropogenic emissions of CO2, geologic storage will be one of the key technologies for achieving this goal. In the standard approach to storage, CO2 is captured from fixed sources such as coal-fired power plants, compressed and injected at supercritical conditions into a suitable target formation. For typical geothermal gradients, suitable formations are found at depths of 800 m (2600 ft) or more, as their temperatures and pressures will be above the critical point of CO2. Trapping the injected CO2 involves one or more mechanisms (Bachu et al. 1994; IPCC 2006): (1) permeability trapping by an impervious confining layer or cap rock; (2) solubility trapping by CO2 dissolution into the aqueous phase in the pore space; (3) mineralogic trapping by chemical reaction of cations with dissolved CO2 to precipitate carbonate minerals; (4) residual phase trapping as the nonwetting CO2 phase becomes disconnected in pores or small clusters of pores; (5) stratigraphic trapping below a formation whose capillary entry pressure is greater than the capillary pressure of the CO2 phase. An intact confining layer is necessary for several trapping mechanisms. However, sedimentary basins often contain geological discontinuities which are potential pathways for leakage through the confining layer. Faults are one such discontinuity and are prevalent in many regions where CO2 storage is likely to be implemented. Wells are a man-made discontinuity, likewise prevalent in likely storage regions. We do not treat them here, but the conceptual model for faults provides a foundation for assessing leakage along wells (Huerta et al. 2008). It is therefore important to examine the consequences if injected CO2 encounters a fault. Figure 1 illustrates the situation of interest. A conductive fault can be a major pathway for the CO2 plume due to its large transfer capacity. However CO2 leaking from the main target formation does not necessarily reach the Earth’s surface. It may not even reach shallower formations of economic interest (mines, hydrocarbon reservoirs, aquifers that serve as underground sources of drinking water (USDW)). Instead, the rising CO2 can be secondarily trapped by shallow subsurface structures, dissolution and residual phase creation (Lindeberg 1997). It can also migrate into permeable formations encountered by the conductive fault. On one hand, this migration attenuates the upward flux. On the other, it spreads the influence of the CO2 across a wider area. The near-surface zone can also attenuate CO2 leaks and decrease CO2 concentration reaching the surface. The attenuation rate is sensitive to the subsurface properties (Oldenburg and Unger 2003). Thus, the effect of a conductive fault on net CO2 storage needs to be analyzed based on the geometric and petrophysical properties of the formation, of the fault and of the overlying permeable layer, and on the boundary conditions (pressure in the storage formation and in the overlying layers). The physics of a CO2 plume rising long vertical distances through the Earth’s crust can be complex (Pruess 2003). Here we present a highly simplified model, motivated by two considerations. One is geological: in many sedimentary basins, a fault is unlikely to be conductive continuously from depth to the shallow subsurface. Leakage to surface or to USDW will involve a sequence of upward (along a fault) and lateral (within a permeable layer) migrations. Thus we will consider moderate conduit lengths of 1000 meters or less. The other consideration is practical: for CO2 storage to be implemented broadly and rapidly enough to mitigate anthropogenic emissions, thousands of storage projects will be needed. Each will have to be permitted by regulators in a streamlined yet robust and transparent way. Unfortunately the physical properties of most storage formations – deep saline aquifers – will be poorly constrained prior to injection. In light of this uncertainty, simple models that allow adequate physics-based risk assessment will be valuable tools for operators, regulators and policymakers. The model presented here was developed to be applicable within the Certification Framework (CF) for geologic storage (Oldenburg and Bryant 2007, Oldenburg et al. 2008). The purpose of the CF is to provide a framework for project proponents, regulators, and the public to analyze the risks of geologic CO2 storage in a simple and transparent way. The risk analysis would be performed to certify the startup and decommissioning of sites for geologic CO2 storage. The CF currently emphasizes risks associated with subsurface processes and excludes compression, transportation, and injection-well leakage risk. The CF is designed to be simple by (1) using proxy concentrations or fluxes, rather than complicated exposure functions, for quantifying impact; (2) using a catalog of pre-computed CO2 injection results (Kumar 2008), and (3) using a simple framework for calculating leakage risk. For transparency, the CF endeavors to be clear and precise in terminology in order to communicate to the full spectrum of stakeholders. One concept of the CF is that leakage occurs along conduits from the storage volume to “compartments” such as hydrocarbon reservoirs or USDW. The risk associated with leakage is the product of the probability of leakage and the impact of that leakage. The flux of CO2 contributes to impact. Thus the goal of the faultleakage model is to estimate flux at an “outlet” of a conductive fault, once CO2 has arrived at the “inlet”. The overarching criteria of the CF and the models within it are simplicity, transparency and acceptability. Here we emphasize the requirement of simplicity. The subsurface data available as input to any numerical model in the CF will always be limited. This limitation is especially troublesome in the case of injection into saline aquifers. Saline aquifers provide large storage capacity but are not well characterized due to the small number of existing wells which could provide geologic information about that particular region of the subsurface. The flow properties of faults are even more uncertain, at least before injection begins. In principle the required properties could be obtained from appropriate measurement campaigns. In practice, measurement will increase the cost of storage, and cost minimization will be a high priority for any greenhouse gas mitigation strategy. Moreover, the most reliable measurements would come from wells drilled into the formation or through the fault. These wells would themselves be potential pathways for leakage. The philosophy of the fault-leakage model is thus to identify the key physical phenomena controlling leakage flux. The sensitivity of the flux to physical parameters provides insight as to which aquifer properties should be measured when designing or monitoring a storage project. In subsequent sections we present a quasi-1D mathematical description of CO2 flux using the “leaky conduit” model. The leaks correspond to permeable formations intersected by the conduit. A “leakoff coefficient” is used to control the rate of leakage. We describe a method for estimating these coefficients from the properties of the formations. To test whether these idealizations are reasonable, we carry out simulations of the full physics of the problem in a 2D domain. Quasi-1D Modeling Approach Assumptions. We assume that the CO2 storage reservoir is located at sufficient depth for the carbon dioxide to be modeled as a slightly compressible fluid. Intersecting this storage volume is a fault, either vertical or at a fixed angle to the storage reservoir (Fig. 1). The most stringent assumption of the model is that we only consider the flow of a single fluid, namely CO2, along the fault. Leakage of CO2 from saline aquifers will involve cocurrent and countercurrent flow of two fluid phases (brine and carbon dioxide), and relative permeability and capillary pressure play critical roles in buoyancy driven flow (e.g. Lakshminarasimhan et al. (2006) and Saadatpoor et al. (2008)). However, for the purposes of quantifying leakage mechanisms, single phase flow of CO2 (Darcy’s law) is assumed to capture the relationship between flux and driving force. We test this idealization later in this paper via two-phase flow simulations carried out using the GEM commercial simulator (Nghiem et al. 2006). We treat the fault as a one-dimensional conduit with spatially varying permeability. Thus the complexities of the fault core and the damage zone surrounding it are averaged into a single array of permeability values. This simplification is extreme but it is consistent with concepts such as the shale gouge ratio (Yielding et al. 1997). Permeability in the direction of the fault throw is likely to be small in regions where sufficient shale is entrained in the fault. The existence of such a region will control CO2 flux along the fault. This is true even in a long fault that elsewhere juxtaposes sand against sand and is very conductive. Accounting for permeability variation along the fault will be important in estimating the risks associated with faults for a geologic storage project. We further assume that the CO2 remains slightly compressible as it rises along the fault. If the CO2 moves slowly enough to equilibrate with the pressure and temperature of the surrounding formations, slight compressibility is a reasonable assumption for two situations. One is for deep leaks, those along faults beginning and ending at depths greater than about 1200 m (3500 ft). The other situation is shallow leaks, along faults less than 600 m (2000 ft) in depth. The density of CO2 for hydrostatic pressure gradient and a representative geothermal gradient (29.2°C/1 km or 1.6°F/100 ft, Tsurface = 15°C or 59°F) is shown as a function of depth in Fig. 2a. The combined variation of pressure and temperature yields an almost linear variation of CO2 density with depth for D > 3500 ft and for D < 2000 ft. The variation in each region can be approximated well with a slightly compressible model with constant compressibility (Table 1). Leaks which rise through intermediate depths, 2000 ft < D < 3500 ft, could also be approximated with a slightly compressible model, and long leakage paths could be modeled by piecing together the appropriate behavior. We do not explore that possibility here. All results are for a single value of compressibility. For simplicity we assume constant viscosity. Figure 2b shows that this is a reasonable assumption if we take a value around 0.000018 Pa-s (0.018 cp) in the “shallow leak” region and 0.00006 Pa-s (0.06 cp) in the “deep leak” region, respectively. This simplification of the phase behavior also assumes that the pressure-temperature profile does not cross the gas/liquid phase boundary. As indicated in Fig. 2c, this assumption is valid for typical geothermal gradients and surface temperatures. The model is “quasi”-1D because CO2 is allowed to enter neighboring strata lateral to the fault (Fig. 3), but flow in these lateral formations is not explicitly modeled. Because we represent the fault as a “leaky conduit”, we refer to flow into these strata as “leakage”. The connotation of “leakage,” that is, flow from the fault vs. flow from the storage formation, should be clear from the context. We will also use the term “attenuation” to indicate CO2 leaving the fault and entering a permeable layer. We model leakage from the fault via a specially designed source term. The lateral formations may be horizontal or inclined. The angle of inclination adds buoyancy to the driving force for leakage into the lateral formation. We assume that the fault contains water initially (implemented by assuming the initial pressure along the fault is hydrostatic). The top of the fault is assumed to be at hydrostatic pressure. The bottom boundary of the fault is set to a pressure above hydrostatic. The value of this over-pressurization is proportional to the height of the column of stored CO2. The pressure in the CO2 phase is taken to be hydrostatic at the base of the column. The thicker the column, the higher the pressure in the CO2 at the location where the fault intersects the storage formation. This high pressure at the bottom of the fault (along with buoyancy) is the driving force for flow along the fault. Mathematical Model. The single-phase flow equation derived from conservation of mass and Darcy’s law is given by Ewing (1983): ( ) ( ) k p g D q t ρ φρ ρ μ ∂ = ∇ ⋅ ∇ − ∇ + ∂ (1) where φ is porosity, ρ is fluid density, p is pore pressure, k is vertical permeability of the fault, μ is fluid viscosity, g is the gravitational constant, D is the depth vector, and q denotes the source or sink term. Mathematically, the boundary and initial conditions described in the previous subsection are modeled as follows. If pz denotes hydrostatic pressure then our initial condition along the fault at initial time t0 is p(z,t0) = pz(z). The pressure at the top of the fault is given by p(ztop, t) = pz(ztop), and the pressure at the bottom of the fault is p(zbottom, t) = pz(zbottom) + X where X is an additional pressure due to the CO2 storage compartment below the fault. Finally the source term is modeled by a vector of “leakoff coefficients” times a potential difference. The coefficients and the potential difference vary along the fault. Physically these leakoff values are only nonzero at depths z corresponding to leaks from the fault into lateral strata. The leakage source term is, therefore, described by ( ) ( ) ( ) ( ) ( ) , , sin leak z q z t q z p z t p z z g ρ γ = − − ⎡ ⎤ ⎣ ⎦ (2) Here in the gravity term, the gradient of the depth vector has been evaluated in the direction of the permeable formation. It is zero for horizontal formations and equal to sin γ for formations inclined at angle γ to the horizontal. Equation 1 is a single equation in two unknowns (density ρ and pressure p). To reduce this equation to one equation in one unknown we make use of the equation of state which relates density to pressure via the differential equation (Chen et al. 2006): f p c t t ρ ρ ∂ ∂ = ∂ ∂ (3) The end result is that our one-dimensional model requires solution of the following equation for fluid pressure along the fault (Chen et al. 2006):